#########################################
################ plots ###############
#########################################

rm(list=ls())
# working directory should be folder extr_weather_rep:
getwd()

# reproducability
print(sessionInfo())
# R version 4.3.0 (2023-04-21 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19045)
# 
# Matrix products: default
# 
# 
# locale:
#   [1] LC_COLLATE=German_Germany.utf8 
# [2] LC_CTYPE=German_Germany.utf8   
# [3] LC_MONETARY=German_Germany.utf8
# [4] LC_NUMERIC=C                   
# [5] LC_TIME=German_Germany.utf8    
# 
# time zone: Europe/Berlin
# tzcode source: internal
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets 
# [6] methods   base     
# 
# other attached packages:
#   [1] rio_1.0.1
# 
# loaded via a namespace (and not attached):
#   [1] vctrs_0.6.3       knitr_1.44        cli_3.6.1        
# [4] gt_0.9.0          xfun_0.40         rematch2_2.1.2   
# [7] rlang_1.1.1       stringi_1.7.12    generics_0.1.3   
# [10] glue_1.6.2        gtExtras_0.5.0    colorspace_2.1-0 
# [13] htmltools_0.5.6   scales_1.2.1      fansi_1.0.4      
# [16] grid_4.3.0        munsell_0.5.0     tibble_3.2.1     
# [19] fontawesome_0.5.2 fastmap_1.1.1     lifecycle_1.0.3  
# [22] compiler_4.3.0    dplyr_1.1.3       pkgconfig_2.0.3  
# [25] rstudioapi_0.15.0 R.oo_1.25.0       R.utils_2.12.2   
# [28] digest_0.6.33     paletteer_1.6.0   R6_2.5.1         
# [31] tidyselect_1.2.0  utf8_1.2.3        pillar_1.9.0     
# [34] magrittr_2.0.3    R.methodsS3_1.8.2 tools_4.3.0      
# [37] gtable_0.3.4      xml2_1.3.5        ggplot2_3.4.3   

# Libraries
library(foreign)
library(tidyverse)
library(rio)
library(gt)
library(gtExtras)
library(raster)
library(sf)
library(ggpubr)
library(ggalluvial)
library(RColorBrewer)
library(ggspatial)
library(ggforce)
library(svglite) # for descr table
library(webshot2)# for descr table


# Anonymise survey data ####
# resp <- import("data/respondents_final.csv")
# #lat: 45.47    46.71    47.15    94.65    47.40 47241.00
# #lng: 5.996   7.447   8.279   8.060   8.731  10.428
# 
# # produce random x and ys for replication data
# resp$lat <- runif(nrow(resp), min = 45.48, max = 47.4)
# resp$lng <- runif(nrow(resp), min = 6, max = 10.4)
# write.csv(resp, "data/respondents_final.csv")
# IDs are also anonymised

# Import anonymised survey data ####
resp <- import("data/respondents_final.csv")

# exclude all respondents that are outside of Switzerland
resp <- dplyr::filter(resp, ID != 39)
resp <- dplyr::filter(resp, ID != 1000)
resp <- dplyr::filter(resp, ID != 1206)
resp <- dplyr::filter(resp, ID != 2117)

## Summary Descriptives of dependent variable ####

# "w4_q15x3_rev", "w1_q16x4_rev"
# "step3_index_w4", "step3_index_w1"
# "w4_q5x9", "w1_q5x9"
# "w4_q29", "w1_q25"
# "w4_q30", "w1_q26"
# "policy_pref_w4", "policy_pref_w1"
# "pid_green_w4", "pid_green_w1"

# code differences

variance <- resp[,FALSE]

variance$diff_step2 <- resp$w4_q15x3_rev - resp$w1_q16x4_rev
variance$diff_step3 <- resp$step3_index_w4 - resp$step3_index_w1
variance$diff_step4_1 <- resp$w4_q5x9 - resp$w1_q5x9
variance$diff_step4_2 <- resp$w4_q29 - resp$w1_q25
variance$diff_step4_3 <- resp$w4_q30 - resp$w1_q26
variance$diff_step4_4 <- resp$policy_pref_w4 - resp$policy_pref_w1

variance$diff_step4_5 <- NA
variance$diff_step4_5[!is.na(resp$pid_green_w4)&!is.na(resp$pid_green_w1)] <-  resp$pid_green_w4[!is.na(resp$pid_green_w4)&!is.na(resp$pid_green_w1)] - 
  resp$pid_green_w1[!is.na(resp$pid_green_w4)&!is.na(resp$pid_green_w1)]


# make table that shows
# A) share of respondents without any variance, missings
# B) hist of changes

# columns
colnames <- c("Variable", "% no change", "n nonmissing", "n change")
variables <- c("Step 2", "Step 3", "Step 4 (MIP)", "Step 4 (Cand. loc.)",
               "Step 4 (Cand. nat.)", "Step 4 (Policy Pref)",
               "Step 4 (green PID)")
nochange <- c(length(variance$diff_step2[variance$diff_step2==0&!is.na(variance$diff_step2==0)]), 
              length(variance$diff_step3[variance$diff_step3==0&!is.na(variance$diff_step3==0)]),
              length(variance$diff_step4_1[variance$diff_step4_1==0&!is.na(variance$diff_step4_1==0)]), 
              length(variance$diff_step4_2[variance$diff_step4_2==0&!is.na(variance$diff_step4_2==0)]),
              length(variance$diff_step4_3[variance$diff_step4_3==0&!is.na(variance$diff_step4_3==0)]), 
              length(variance$diff_step4_4[variance$diff_step4_4==0&!is.na(variance$diff_step4_4==0)]),
              length(variance$diff_step4_5[variance$diff_step4_5==0&!is.na(variance$diff_step4_5==0)]))
nonmissing <- c(length(variance$diff_step2[!is.na(variance$diff_step2)]), 
              length(variance$diff_step3[!is.na(variance$diff_step3)]),
              length(variance$diff_step4_1[!is.na(variance$diff_step4_1)]), 
              length(variance$diff_step4_2[!is.na(variance$diff_step4_2)]),
              length(variance$diff_step4_3[!is.na(variance$diff_step4_3)]), 
              length(variance$diff_step4_4[!is.na(variance$diff_step4_4)]),
              length(variance$diff_step4_5[!is.na(variance$diff_step4_5)]))
med_var <- c(median(variance$diff_step2[variance$diff_step2!=0], na.rm=T),
             median(variance$diff_step3[variance$diff_step3!=0], na.rm=T),
             median(variance$diff_step4_1[variance$diff_step4_1!=0], na.rm=T),
             median(variance$diff_step4_2[variance$diff_step4_2!=0], na.rm=T),
             median(variance$diff_step4_3[variance$diff_step4_3!=0], na.rm=T),
             median(variance$diff_step4_4[variance$diff_step4_4!=0], na.rm=T),
             median(variance$diff_step4_5[variance$diff_step4_5!=0], na.rm=T))
mean_var <- c(mean(variance$diff_step2[variance$diff_step2!=0], na.rm=T),
             mean(variance$diff_step3[variance$diff_step3!=0], na.rm=T),
             mean(variance$diff_step4_1[variance$diff_step4_1!=0], na.rm=T),
             mean(variance$diff_step4_2[variance$diff_step4_2!=0], na.rm=T),
             mean(variance$diff_step4_3[variance$diff_step4_3!=0], na.rm=T),
             mean(variance$diff_step4_4[variance$diff_step4_4!=0], na.rm=T),
             mean(variance$diff_step4_5[variance$diff_step4_5!=0], na.rm=T))

             
nchange <- nonmissing-nochange 
nochange_perc <- nochange/(nochange+nchange) * 100
change_perc <- nchange/(nochange+nchange) *100
min <- apply(variance, 2, min, na.rm =T)
max <- apply(variance, 2, max, na.rm =T)

variance_sum <- cbind(variables, min, max, nonmissing, nochange_perc,  change_perc, med_var, mean_var) %>% data.frame()
variance_sum$nonmissing <- as.numeric(variance_sum$nonmissing)
variance_sum$nochange_perc <- as.numeric(variance_sum$nochange_perc )
variance_sum$change_perc <- as.numeric(variance_sum$change_perc )
variance_sum$min <- as.numeric(variance_sum$min)
variance_sum$max <- as.numeric(variance_sum$max)
variance_sum$med_var <- as.numeric(variance_sum$med_var)
variance_sum$mean_var <- as.numeric(variance_sum$mean_var)
variance_sum$dens_data <- apply(variance, 2, list)
# list has one layer too many.
variance_sum$dens_data2[1] <- variance_sum$dens_data[[1]]
variance_sum$dens_data2[2] <- variance_sum$dens_data[[2]]
variance_sum$dens_data2[3] <- variance_sum$dens_data[[3]]
variance_sum$dens_data2[4] <- variance_sum$dens_data[[4]]
variance_sum$dens_data2[5] <- variance_sum$dens_data[[5]]
variance_sum$dens_data2[6] <- variance_sum$dens_data[[6]]
variance_sum$dens_data2[7] <- variance_sum$dens_data[[7]]
variance_sum$dens_data <- variance_sum$dens_data2
variance_sum$dens_data2 <- NULL

table_desc <- variance_sum %>% 
  gt() %>%
  # histogram
  gtExtras::gt_plt_dist(
  dens_data,
  type = "histogram",
  line_color = "#474747FF",
  fill_color = "#474747FF",
  bw = 0.75,
  same_limit = TRUE
  ) %>% 
  # format decimals
  fmt_number(columns = min:max, decimals = 0) %>%
  fmt_number(columns = nochange_perc:change_perc, decimals = 2) %>%
  fmt_number(columns = med_var:mean_var, decimals = 2) %>%
  # header
tab_header(
  title = md("**Summary of Dependent Variables**"),
  subtitle = md("Based on first differences for each step,
                including information on level of variation within individuals")) %>%
  tab_spanner(label = "Summary statistics",
    columns = min:nonmissing) %>%
  tab_spanner(label = "Graphics",
    columns = dens_data)  %>%
  tab_spanner(label = "Variation",
    columns = nochange_perc:change_perc) %>%
  tab_spanner(label = "Intensity of change",
              columns = med_var:mean_var) %>%
  # change column names to appear in the table
  cols_label(variables = "Variables within respondents",
             nochange_perc = html("no (%)"),
             change_perc = html("yes (%)"),
             med_var = html("median"),
             mean_var = html("mean"),
             dens_data = "Histogram") %>%
  # change some defaults on table appearance
  tab_options(data_row.padding = px(6),
    heading.padding = px(0),
    column_labels.padding = px(0),
    footnotes.padding = px(0),
    table_body.hlines.width = px(0.5),
    table_body.hlines.color = "#474747FF",
    table_body.vlines.width = px(0.5),
    table_body.vlines.color = "#474747FF",
    row_group.border.bottom.width = px(0.5),
    row_group.border.bottom.color = "#474747FF",
    row_group.border.top.width = px(0.5),
    row_group.border.top.color = "#474747FF",
    column_labels.border.bottom.width = px(0.5),
    column_labels.border.bottom.color = "#474747FF",
    column_labels.border.top.width = px(0.5),
    column_labels.border.top.color = "#474747FF",
    table_body.border.bottom.width = px(0.5),
    table_body.border.bottom.color = "#474747FF",
    table_body.border.top.width = px(0.5),
    table_body.border.top.color = "#474747FF",
    table.border.bottom.width = px(1.5),
    table.border.bottom.color = "#474747FF",
    table.border.top.width = px(1.5),
    table.border.top.color = "#474747FF",
    heading.border.bottom.width = px(1.5),
    heading.border.bottom.color = "#474747FF",
    stub.border.width = px(1.5),
    stub.border.color = "#474747FF",
    table.background.color = "white",
    table.font.color = "#474747FF",
    #table.font.names = "Open Sans",
    table.font.weight = "normal",
    table.font.size = "14px",
    heading.subtitle.font.size = "14px",
    column_labels.font.weight = "bold",
    footnotes.font.size = "12px") %>%
  # set alignment as per wish
  cols_align(align = "left",
    columns = everything()) %>%
  opt_align_table_header(align = "left") %>%
  # set column widths
  cols_width(variables:mean_var ~ px(150))%>%
  tab_footnote(footnote = "calculated only for those cases with non-zero variation",
  locations = cells_column_spanners(spanners = "Intensity of change"))
  
gtsave(table_desc, 
  "plots/table_dv.png", expand = 10)

rm(variance); rm(variance_sum); rm(change_perc); rm(max)
rm(mean_var); rm(med_var); rm(min); rm(nchange); rm(nochange)
rm(nochange_perc); rm(nonmissing); rm(variables); rm(colnames)
rm(table_desc)

## Summary Descriptives of independent variable ####

# get labels of weather variables
operat <- import("data/operationalisations.xlsx")
ind_var_labels <- operat[,2] %>% as.vector()
rm(operat)
n_start <- which(names(resp) == "temp_extreme_w1w4_percent")
n_end <- which(names(resp) == "dis_binary_w1w4_summer")
ind_var <- names(resp[,c(n_start:n_end)])

# create empty vectors to fill
nonmissing <- rep(NA, length(ind_var))

# fill with values
for (i in 1:length(ind_var)) {
  variable <- which(names(resp)==ind_var[i])
  pick_variable <- resp[, variable]
  nonmissing[i] <- length(pick_variable[!is.na(pick_variable)])
  }

min <- apply(resp[,n_start:n_end], 2, min, na.rm =T)
quart1 <- apply(resp[,n_start:n_end], 2, quantile, 0.25)
median <- apply(resp[,n_start:n_end], 2, median, na.rm =T)
mean <- apply(resp[,n_start:n_end], 2, mean, na.rm =T)
quart3 <- apply(resp[,n_start:n_end], 2, quantile, 0.75)
max <- apply(resp[,n_start:n_end], 2, max, na.rm =T)

variance_sum <- cbind(1:length(ind_var_labels),ind_var_labels, min, quart1, median, mean, quart3, max, 
                      nonmissing) %>% data.frame()
variance_sum$nonmissing <- as.numeric(variance_sum$nonmissing)
variance_sum$min <- as.numeric(variance_sum$min)
variance_sum$quart1 <- as.numeric(variance_sum$quart1)
variance_sum$median <- as.numeric(variance_sum$median)
variance_sum$mean <- as.numeric(variance_sum$mean)
variance_sum$quart3 <- as.numeric(variance_sum$quart3)
variance_sum$max <- as.numeric(variance_sum$max)

resp$dis_binary_w1w4 <- as.numeric(resp$dis_binary_w1w4)
resp$dis_binary_w1w4_summer <- as.numeric(resp$dis_binary_w1w4_summer)

variance_sum$dens_data <- apply(resp[,n_start:n_end], 2, list)
# list has one layer too many.
for (i in 1:nrow(variance_sum)){
  variance_sum$dens_data2[i] <- variance_sum$dens_data[[i]]
}
variance_sum$dens_data <- variance_sum$dens_data2
variance_sum$dens_data2 <- NULL

table_ind_var <- variance_sum %>% 
  gt() %>%
  gtExtras::gt_plt_dist(
    dens_data,
    type = "histogram",
    line_color = "#474747FF",
    fill_color = "#474747FF",
    bw = 0.75,
    same_limit = TRUE,
    trim = F
  ) %>%
  # format decimals
  fmt_number(columns = min:max, decimals = 2) %>%
  # header
  tab_header(
    title = md("**Summary of Weather Variables**"),
    subtitle = md("Summary statistics")) %>%
  # tab_spanner(label = "Summary statistics",
  #             columns = min:nonmissing) %>%
  # tab_spanner(label = "Graphics",
  #             columns = dens_data)  %>%
  # # change column names to appear in the table
  cols_label(V1 = "No.",
             ind_var_labels = "Weather variable",
             quart1 = "25%",
             quart3 = "75%",
             dens_data = "Histogram",
             nonmissing = "n") %>%
  # change some defaults on table appearance
  tab_options(data_row.padding = px(0),
              heading.padding = px(0),
              column_labels.padding = px(0),
              footnotes.padding = px(0),
              table_body.hlines.width = px(0.5),
              table_body.hlines.color = "#474747FF",
              table_body.vlines.width = px(0.5),
              table_body.vlines.color = "#474747FF",
              row_group.border.bottom.width = px(0.5),
              row_group.border.bottom.color = "#474747FF",
              row_group.border.top.width = px(0.5),
              row_group.border.top.color = "#474747FF",
              column_labels.border.bottom.width = px(0.5),
              column_labels.border.bottom.color = "#474747FF",
              column_labels.border.top.width = px(0.5),
              column_labels.border.top.color = "#474747FF",
              table_body.border.bottom.width = px(0.3),
              table_body.border.bottom.color = "#474747FF",
              table_body.border.top.width = px(0.3),
              table_body.border.top.color = "#474747FF",
              table.border.bottom.width = px(1.5),
              table.border.bottom.color = "#474747FF",
              table.border.top.width = px(1.5),
              table.border.top.color = "#474747FF",
              heading.border.bottom.width = px(1.5),
              heading.border.bottom.color = "#474747FF",
              stub.border.width = px(1.5),
              stub.border.color = "#474747FF",
              table.background.color = "white",
              table.font.color = "#474747FF",
              #table.font.names = "Open Sans",
              table.font.weight = "normal",
              table.font.size = "14px",
              heading.subtitle.font.size = "14px",
              column_labels.font.weight = "bold",
              footnotes.font.size = "12px") %>%
  # set alignment as per wish
  cols_align(align = "left",
             columns = everything()) %>%
  opt_align_table_header(align = "left") %>%
  # set column widths
  cols_width(min:nonmissing ~ px(200)) %>%
  cols_width(ind_var_labels ~ px(470))
  
gtsave(table_ind_var, 
       "plots/table_ind_var.png", 
       expand = 15)

rm(table_ind_var); rm(variance_sum); rm(i); rm(ind_var); rm(ind_var_labels)
rm(max); rm(mean); rm(median); rm(min); rm(n_end); rm(n_start)
rm(nonmissing); rm(pick_variable); rm(quart1); rm(quart3); rm(variable)

# Overview of Operationalisations in table

operat <- import("data/operationalisations.xlsx")
operat[is.na(operat)] <- "" 
operat_tab1 <- operat[1:18,] %>% gt() %>%
  # set alignment as per wish
  cols_align(align = "left",
             columns = everything()) %>%
  opt_align_table_header(align = "left") %>%
  # set column widths
  cols_width(starts_with("Descr") ~ px(550)) %>%
  cols_width(starts_with("Variable") ~ px(150))
gtsave(operat_tab1, "plots/operat_p1.png", 
       expand = 10)

operat_tab2 <- operat[19:34,] %>% gt() %>%
  # set alignment as per wish
  cols_align(align = "left",
             columns = everything()) %>%
  opt_align_table_header(align = "left") %>%
  # set column widths
  cols_width(starts_with("Descr") ~ px(550)) %>%
  cols_width(starts_with("Variable") ~ px(150))
gtsave(operat_tab2, "plots/operat_p2.png", 
       expand = 10)

rm(operat); rm(operat_tab1); rm(operat_tab2)

## Stability of causal chain over steps ####

stab <- import("data/chains_of_steps.xlsx")
stab_tab1 <- stab[,1:5] %>% gt() %>%
  # set alignment as per wish
  cols_align(align = "left",
             columns = everything()) %>%
  opt_align_table_header(align = "left") %>%
  # set column widths
  cols_width(everything() ~ px(180)) %>%
  cols_width(starts_with("Step") ~ px(250)) 

gtsave(stab_tab1, "plots/stab_tab1.png", 
       expand = 10)

rm(stab); rm(stab_tab1)

# Geospatial plots ####

# Shapefile municipalities
# download this data from here: https://www.swisstopo.admin.ch/de/landschaftsmodell-swissboundaries3d
# this analysis is based on 2020 data
# save it to data/municipalities
muni <- shapefile("data/municipalities/swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET.shp")
# only keep information on muni number
muni <- muni[,"BFS_NUMMER"]

muni$ID <- 1:nrow(muni)
# there are 2350 municipalities

extent(muni)
# class      : Extent 
# xmin       : 2485410 
# xmax       : 2833842 
# ymin       : 1075268 
# ymax       : 1295934
# correct extent: https://epsg.io/2056
muni$srph_municipal <- muni$BFS_NUMMER
muni$BFS_NUMMER <- NULL
muni$ID <- NULL
crs(muni) <- "+init=epsg:2056"

# boundaries Switzerland for background
# get boundaries from here: https://gadm.org/download_country.html (level 0 Switzerland)
# save in data folder
CHE <- read_rds("data/gadm36_CHE_0_sp.rds")
CHE <- st_as_sf(CHE)
CHE <- st_transform(CHE, "+init=epsg:2056")

# lakes 
# download this data from here: https://www.swisstopo.admin.ch/de/hoehenmodell-swissbathy3d#Technische-Details-pro-See
# this analysis is based on 2019 data that is no longer public, but should also work with new data
# save it to data/lakes
lakes <- read_sf("data/lakes/g2s19.shp")
lakes <- lakes[,2]
lakes_sf <- st_transform(lakes, "+init=epsg:2056")
# cut to borders switzerland
lakes_sf <- st_intersection(lakes_sf, CHE)[,1]
rm(lakes)

# rivers
# download (very similar, but not identical data) here:
# https://www.geocat.ch/geonetwork/srv/ger/catalog.search#/metadata/30b22e49-7ef5-4fa6-b957-7b38cc3beb27
# save it to data/rivers
rivers1 <- read_sf("data/rivers/k4flusyyyymmdd55_ch2007.shp")
rivers2 <- read_sf("data/rivers/k4flusyyyymmdd11_ch2007.shp")
rivers3 <- read_sf("data/rivers/k4flusyyyymmdd22_ch2007.shp")
rivers4 <- read_sf("data/rivers/k4flusyyyymmdd33_ch2007.shp")
rivers5 <- read_sf("data/rivers/k4flusyyyymmdd44_ch2007.shp")
# combine to 1
rivers <- rbind(rivers1, rivers2, rivers3, rivers4, rivers5)[,1]
rm(rivers1); rm(rivers2); rm(rivers3); rm(rivers4); rm(rivers5); 
# cut to borders switzerland
rivers_sf <- st_intersection(rivers, CHE)[,1]
rm(rivers)

# productive cantons
# productive cantons no longer available for download,
# but maps of settlement very similar: https://data.geo.admin.ch/browser/index.html#/collections/ch.are.siedlung?.language=en&.asset=asset-siedlung.zip
# download and save to data/productive_cantons
canton_prod_geo <- read_sf("data/productive_cantons/K4kant19970101vf_ch2007Poly.shp")
# plot(canton_prod_geo[,1])

# Distribution of all treatment variables ####

# disasters: 6 operationalisations
ggplot(resp, aes(dis_extreme_w1w4_percent)) + geom_density() +
  geom_density(aes(dis_extreme_w1w4_summer_percent)) +
  geom_density(aes(dis_extreme_w1w4_25y_percent)) +
  geom_density(aes(dis_extreme_w1w4_25y_summer_percent)) +
  geom_density(aes(dis_binary_w1w4)) +
  geom_density(aes(dis_binary_w1w4_summer)) +
  theme_bw() + labs(x = "Extremeness of events taken place", y = "density", 
                    title = "Distribution of extreme events between the two survey waves")

ggsave("plots/distr_events.png", 
       width = 12, height = 6, dpi = 400)

# rain: 8 operationalisations
ggplot(resp, aes(rain_extreme_w1w4_percent)) + geom_density() +
  geom_density(aes(rain_extreme_w1w4_summer_percent)) +
  geom_density(aes(rain_extreme_w1w4_25y_percent)) +
  geom_density(aes(rain_extreme_w1w4_25y_summer_percent)) +
  geom_density(aes(rain_extreme_w1w4_drought_percent)) +
  geom_density(aes(rain_extreme_w1w4_drought_summer_percent)) +
  geom_density(aes(rain_extreme_w1w4_spi_percent)) +
  geom_density(aes(rain_extreme_w1w4_spi_summer_percent)) +
  theme_bw() + scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(limits = c(0,0.6), expand = c(0, 0)) +
  labs(x = "percent of extreme days", y = "density", 
       title = "Distribution of extreme rain between the two survey waves",
       caption = "(x-axis cropped at 60)")
ggsave("plots/distr_rain.png", 
       width = 12, height = 6, dpi = 400)


# temp relative: 12 operationalisations

temp_rel <- ggplot(resp, aes(temp_extreme_w1w4_summer_percent)) + geom_density() +
  geom_density(aes(temp_extreme_w1w4_percent)) +
  geom_density(aes(temp_extreme_w1w4_25y_percent)) +
  geom_density(aes(temp_extreme_w1w4_25y_summer_percent)) +
  geom_density(aes(temp_mean_extreme_w1w4_percent)) +
  geom_density(aes(temp_mean_extreme_w1w4_summer_percent)) +
  geom_density(aes(temp_mean_extreme_w1w4_25y_percent)) +
  geom_density(aes(temp_mean_extreme_w1w4_25y_summer_percent)) +
  geom_density(aes(temp_extreme_90_w1w4_percent)) +
  geom_density(aes(temp_extreme_90_w1w4_summer_percent)) +
  geom_density(aes(temp_extreme_90_w1w4_25y_percent, fill= "red", alpha=0.3)) +
  geom_density(aes(temp_extreme_90_w1w4_25y_summer_percent)) +
  theme_bw() + scale_x_continuous(expand = c(0, 0),
                                  limits = c(0, 40)) + 
  scale_y_continuous(limits = c(0,0.61), expand = c(0, 0)) +
  labs(x = "percent of extreme days", y = "density", 
       title = "Distribution of relative measures of extreme temperatures",
       caption = "(x-axis cropped at 40)") +
  annotate("text", x = 22, y = 0.33, label = "IPCC",
           color = "#330000") +
  geom_segment(aes(x = 22, y = 0.31, xend = 20, yend = 0.25),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "#330000") +
  theme(legend.position = "none")

# temp absolute: 8 operationalisations

temp_abs <- ggplot(resp, aes(temp_extreme_w1w4_heat25_percent)) + geom_density() +
  geom_density(aes(temp_extreme_w1w4_heat25_summer_percent)) +
  geom_density(aes(temp_extreme_w1w4_heat27_percent, fill= "red", alpha=0.3)) +
  geom_density(aes(temp_extreme_w1w4_heat27_summer_percent)) +
  geom_density(aes(temp_extreme_w1w4_heat30_percent)) +
  geom_density(aes(temp_extreme_w1w4_heat30_summer_percent)) +
  geom_density(aes(temp_extreme_w1w4_heat35_percent)) +
  geom_density(aes(temp_extreme_w1w4_heat35_summer_percent)) +
  theme_bw() + scale_x_continuous(expand = c(0, 0),
                                  limits = c(0, 20)) + 
  scale_y_continuous(limits = c(0,0.6), expand = c(0, 0)) +
  annotate("text", x = 4.5, y = 0.58, label = "MeteoSwiss",
           color = "#330000") +
  geom_segment(aes(x = 2, y = 0.565, xend = 1, yend = 0.50),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "#330000") + 
  labs(x = "percent of extreme days", y = "density", 
       title = "Distribution of absolute measures of extreme temperatures",
       caption = "(x-axis cropped at 20)")+
  theme(legend.position = "none")

# save combined plot distributions
ggarrange(temp_rel, NULL, temp_abs, nrow = 1, 
          widths = c(1, 0, 1))
ggsave("plots/distr_temp.png", 
       width = 12, height = 5, dpi = 400)

rm(temp_rel); rm(temp_abs)

# Dependent Variables ####

# Climate Belief (Step 2)
table(resp$w1_q16x4_rev, resp$w4_q15x3_rev)

w1 <- c(1:5, 1:5, 1:5, 1:5, 1:5)
w4 <- rep(1:5, each = 5)
values <- c(table(resp$w1_q16x4_rev, resp$w4_q15x3_rev))
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

# check whether right format for alluvial plot
is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  ggtitle("Step 2 (Atribution) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/climate_belief_w1w4.png",
       width = 6, height = 4, dpi = 400)

# description:
# mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$w1_q16x4_rev)
summary(resp$w4_q15x3_rev)
# change positive vs negative
resp$change <- resp$w4_q15x3_rev - resp$w1_q16x4_rev  
table(resp$change > 0)
table(resp$change < 0)
table(resp$change == 0)

# Climate Impacts
table(round(resp$step3_index_w1,0), 
      round(resp$step3_index_w4,0))

w1 <- rep(c(1:5), 5)

w4 <- rep(c(1:5), each = 5)
values <- c(table(round(resp$step3_index_w1, 0), 
                  round(resp$step3_index_w4, 0)))
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  ggtitle("Step 3 (Emotional relation, willingness to act) over time")+
  labs(caption = "rounded to full numbers") + 
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/climate_impact_w1w4.png",
       width = 6, height = 5, dpi = 400)

# description:
#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$step3_index_w1)
summary(resp$step3_index_w4)
# change positive vs negative
resp$change <- resp$step3_index_w4 - resp$step3_index_w1
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# MIP over time
table(resp$w1_q5x9, resp$w4_q5x9)
w1 <- c(0, 0, 1, 1)
w4 <- c(0:1, 0:1)
values <- c(1044, 745, 247, 866)
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  scale_fill_brewer(type = "qual", palette = 6) +
  ggtitle("Step 4 (MIP) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/MIP_w1w4.png",
       width = 5, height = 3, dpi = 400)

# description:
#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$w1_q5x9)
summary(resp$w4_q5x9)
# change positive vs negative
resp$change <- resp$w4_q5x9 - resp$w1_q5x9
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# Green local candidate
table(resp$w1_q25, resp$w4_q29)
w1 <- c(1:5, 1:5, 1:5, 1:5, 1:5)
w4 <- rep(c(1:5), each = 5)
values <- c(table(resp$w1_q25, resp$w4_q29))
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  ggtitle("Step 4 (local green candidate) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/green_cand_loc_w1w4.png",
       width = 6, height = 4, dpi = 400)

# description:
#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$w1_q25)
summary(resp$w4_q29)
# change positive vs negative
resp$change <- resp$w4_q29 - resp$w1_q25
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# Green national candidate
table(resp$w1_q26, resp$w4_q30)
w1 <- c(1:5, 1:5, 1:5, 1:5, 1:5)
w4 <- rep(c(1:5), each = 5)
values <- c(table(resp$w1_q26, resp$w4_q30))
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  ggtitle("Step 4 (national green candidate) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/green_cand_nat_w1w4.png",
       width = 6, height = 4, dpi = 400)

# description:
#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$w1_q26)
summary(resp$w4_q30)
# change positive vs negative
resp$change <- resp$w4_q30 - resp$w1_q26
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# Policy Preferences
table(round(resp$policy_pref_w1,0), round(resp$policy_pref_w4, 0))
w1 <- c(1:5, 1:5, 1:5, 1:5, 1:5)
w4 <- rep(c(1:5), each = 5)
values <- c(table(round(resp$policy_pref_w1,0), round(resp$policy_pref_w4, 0)))
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  ggtitle("Step 4 (policy preferences) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/pol_pref_w1w4.png",
       width = 6, height = 4, dpi = 400)

#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$policy_pref_w1)
summary(resp$policy_pref_w4)
# change positive vs negative
resp$change <- resp$policy_pref_w4 - resp$policy_pref_w1
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# Green PID
table(resp$pid_green_w1, resp$pid_green_w4)
w1 <- c(0, 0, 1, 1)
w4 <- c(0:1, 0:1)
values <- c(1243, 172, 53, 310)
resp_climate_wide <- base::as.data.frame(cbind(w1, w4, values), row.names=F)

is_alluvia_form(resp_climate_wide, axes = 1:3, silent = TRUE)

ggplot(resp_climate_wide,
       aes(y = values, axis1 = w1, axis2 = w4)) +
  geom_alluvium(aes(fill = as.factor(values)), width = 1/12) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("w1", "w4"), expand = c(.05, .05)) +
  scale_fill_brewer(type = "qual", palette = 6) +
  ggtitle("Step 4 (Green PID) over time") +
  guides(fill=guide_legend(title="n per group")) +
  theme_void()

ggsave("plots/green_PID_w1w4.png",
       width = 5, height = 3, dpi = 400)

#mean 2018, mean 2019; share of respondents "greening"&reverse
summary(resp$pid_green_w1)
summary(resp$pid_green_w4)
# change positive vs negative
resp$change <- resp$pid_green_w4 - resp$pid_green_w1
table(resp$change > 0)
table(resp$change < 0)
table(resp$change==0)

# Descriptive weather plots (Figure 4) #####

# turn into shapefile
# numeric, exclude cases with missing coordinates
resp$lat <- as.numeric(resp$lat)
resp$lng <- as.numeric(resp$lng)
resp <- resp[!is.na(resp$lat),] # exclude 0 cases
resp <- resp[resp$lat > 35 & resp$lat < 40000,]

resp_spatial <- SpatialPointsDataFrame(coords = resp[c("lng", "lat")],
                                       data = resp,
                                       proj4string = CRS("+init=epsg:4326"))

# turn into NEW swiss crs
resp_spatial <- spTransform(resp_spatial, CRS("+init=epsg:2056"))

# check whether extent matches municipalities:
extent(resp_spatial) 
extent(muni)

# turn into simple feature for ggplot 
resp_sf <- st_as_sf(resp_spatial)
resp_sf <- st_transform(resp_sf, "+init=epsg:2056")

# exclude cases that are not in switzerland
resp_sf_in_CHE <- st_contains(CHE, resp_sf)
resp_sf_in_CHE <- as.data.frame(resp_sf_in_CHE)
resp_sf$in_CHE <- 0
resp_sf$in_CHE[resp_sf_in_CHE[,2]] <- 1
resp_sf <- resp_sf[resp_sf$in_CHE==1,]
resp_sf$in_CHE <- NULL; rm(resp_sf_in_CHE)
# 0 cases excluded
# 2933 cases remaining

# percent of extreme days between w1 and w4: 10y, max 90 #
# display.brewer.all()
nb.cols <- length(unique(as.factor(round(resp_sf$temp_extreme_90_w1w4_25y_percent,0))))
mycolors <- colorRampPalette(rev(brewer.pal(10, "Spectral")))(nb.cols)

mybreaks <- c(summary(round(resp_sf$temp_extreme_90_w1w4_25y_percent,0))[c(1,2,3,5,6)]) %>% as.character()

# Plot percentage of extreme temp days between w1w4
temp_90_25 <- ggplot() +
  # backgroundmap
  geom_sf(data = CHE, color="darkgrey", fill = "grey70") + 
  # background: productive cantons
  geom_sf(data = canton_prod_geo, color="darkgrey", fill = "lightgrey") + 
  # lakes and rivers
  geom_sf(data = lakes_sf, color="darkgrey", fill = "steel blue 3") + 
  geom_sf(data = rivers_sf, color ="steel blue 3") + 
  # disaster data
  geom_sf(data = resp_sf,  
          size = 1, aes(col = as.factor(round(temp_extreme_90_w1w4_25y_percent,0)))) +
  labs(colour = "Percent") +
  scale_color_manual(values = mycolors) + #,
                     #breaks = mybreaks) + 
  ggtitle("       Percent of extreme temperature days", subtitle = "         max temp > 90%, 25 year reference") +
  annotation_scale(location = "bl", width_hint = 0.5) + 
  theme_void() +
  geom_circle(aes(x0 = 2715422, y0 = 1100937, r = 30000),
              inherit.aes = FALSE, color = "#330000") +
  annotate("text", x = 2780422, y = 1080937, label = "Ticino", color = "#330000") +
  geom_curve(
    aes(x = 2780422, y = 1085937, xend = 2750422, yend = 1100937),
    arrow = arrow(length = unit(0.03, "npc")),
    color = "#330000"
  )


# percentage of days over 27 degrees
# display.brewer.all()
nb.cols <- length(unique(as.factor(round(resp_sf$temp_extreme_w1w4_heat27_percent,0))))
mycolors <- colorRampPalette(rev(brewer.pal(10, "Spectral")))(nb.cols)
mybreaks <- c(summary(round(resp_sf$temp_extreme_w1w4_heat27_percent,0))[c(1,2,3,5,6)]) %>% as.character()

# Plot percentage of extreme temp days between w1w4
temp_27 <- ggplot() +
  # backgroundmap
  geom_sf(data = CHE, color="darkgrey", fill = "grey70") + 
  # background: productive cantons
  geom_sf(data = canton_prod_geo, color="darkgrey", fill = "lightgrey") + 
  # lakes and rivers
  geom_sf(data = lakes_sf, color="darkgrey", fill = "steel blue 3") + 
  geom_sf(data = rivers_sf, color ="steel blue 3") + 
  # disaster data
  geom_sf(data = resp_sf,  
          size = 1, aes(col = as.factor(round(temp_extreme_w1w4_heat27_percent,0)))) +
  labs(colour = "Percent") +
  scale_color_manual(values = mycolors) + #,
                     #breaks = mybreaks) + 
  ggtitle("       Percent of extreme temperature days", subtitle = "         mean temp >= 27 degrees") +
  annotation_scale(location = "bl", width_hint = 0.5) + 
  theme_void() +
  annotate("text", x = 2570422, y = 1285937, label = "Basel",
           color = "#330000") +
  geom_curve(
    aes(x = 2585422, y = 1285937, xend = 2610422, yend = 1275937),
    arrow = arrow(length = unit(0.03, "npc")),
    curvature = -0.5,
    color = "#330000")

# save combined plot
ggarrange(temp_90_25, NULL, temp_27, nrow = 1, 
          widths = c(1, 0, 1))
ggsave("plots/temp_90_25+temp_27.png", 
       width = 12, height = 4, dpi = 400)